ApJ, accepted 

Preprint typeset using I^T^^ style cmulateapj v. 12/14/05 



A FAST NEW PUBLIC CODE FOR COMPUTING PHOTON ORBITS IN A KERR SPACETIME 

Jason Dexter 

Department of Physics, University of Washington, Seattle, WA 98195-1560, USA 



o 
o 

(N 



o 



> 
o 

p 

rn 
O 

ON 

o 



X 



Eric Agol 

Department of Astronomy, University of Washington, Box 351580, Seattle, WA 98195, USA 

ApJ, accepted 

ABSTRACT 

Relativistic radiative transfer problems require the calculation of photon trajectories in curved 
spacetime. We present a novel technique for rapid and accurate calculation of null geodesies in the 
Kerr metric. The equations of motion from the Hamilton-Jacobi equation are reduced directly to 
Carlson's elliptic integrals, simplifying algebraic manipulations and allowing all coordinates to be 
computed semi-analytically for the first time. We discuss the method, its implementation in a freely 
available FORTRAN code, and its application to toy problems from the literature. 
Subject headings: accretion — black hole physics — radiative transfer — relativity 



1. INTRODUCTION 

Efficient and accurate computation of null geodesies 
in the vicinity of spinning black holes is important for 
studies of active galaxies, X-ray binaries, and other ac- 
creting black hole systems. The radiated flux from ac- 
cretion disks mostly originates in the innermost radii, 
where relativistic effects are important for understand- 
ing observations. Proper calculation of the bending of 
light requires integration along rays (Broderick 2006). In 
general, propagation through the plasma will influence 
the photon trajectories, leading to non-geodesic paths 
(Broderick & Blandford 2003, 2004). However, these ef- 
fects are mostly important at low frequencies, compa- 
rable to the expected plasma and cyclotron frequency. 
When plasma effects can be neglected, the rays corre- 
spond to null geodesies, and these circumstances are as- 
sumed throughout this paper. 

The flrst applications of general relativistic radiative 
transfer to accreting systems were of two main types. 
Cunningham (1975) packaged all radiative effects for op- 
tically thick, geometrically thin disks as a transfer func- 
tion to go from local emissivity to that observed at in- 
finity. Luminet (1979) used the simple relationships be- 
tween impact parameters at infinity and constants of the 
motion to shoot rays backwards in time from an ob- 
server's photographic plate to the object under study. 
More recently, Viergutz (1993) and Beckwith & Done 
(2005) considered the so-called emitter-observer prob- 
lem. That is, given locations of the emitter and the 
observer, determine the constants of the motion for null 
geodesies connecting the two. This approach is much 
more efficient when the source is highly localized, such as 
an orbiting star or hotspot. Here, backwards ray shoot- 
ing is impractical since most of the rays miss the target. 

Such techniques have been applied to the study of emis- 
sion lines and spectra from active galactic nuclei (AGN) 
accretion disks and tori (Cadez et al. 1998; Wu & Wang 
2007) as well as their quasi-periodic oscillations (QPOs) 
(Schnittman et al. 2006). Li et al. (2005) used a ray 
tracing approach to study the spectra of X-ray binaries. 
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Noble et al. (2007) created images of galactic center black 
hole candidate Sagittarius A* (Sgr A*) using axisymmet- 
ric general relativistic MHD (GRMHD) simulations, and 
Bromley et al. (2001) studied its polarization from a sim- 
plified accretion model. Broderick & Loeb (2006) mod- 
eled the frequency dependence of its centroid position, 
and Reid et al. (2008) used ray tracing to compare hot 
spot accretion models with the observed astrometric mo- 
tion of its mean position as a function of wavelength. Fi- 
nally, although the spacetime surrounding neutron stars 
only asymptotically approaches the Kerr metric, using 
its null geodesies for ray tracing has still found applica- 
tion in modeling spectra of neutron stars (Braje et al. 
2000). 

Despite all of this work, numerical integration of Kerr 
null geodesies is computationally expensive in certain ap- 
phcations. Ranch & Blandford (1994) (hereafter RB94) 
described a method for calculating null geodesies in the 
Kerr metric semi-analytically using the Hamilton-Jacobi 
formulation of the equations of motion and used it to 
study the primary caustic. Bozza (2008) used a simi- 
lar method to investigate caustics of all orders, build- 
ing on earlier analytic work (Bozza 2002). Fanton et al. 
(1997) used a fast analytic version for creating line pro- 
files and accretion disk images, and Agol (1997) applied 
this method to the case of polarization from thin disk 
accretion. Falcke et al. (2000) went on to use this code 
along with a simple model for the Galactic center black 
hole to create images of its accretion fiow. 

All of this work used Legendre's formulation of elliptic 
integrals (e.g., Abramowitz & Stegun 1965), and treated 
the and t coordinates numerically, if at all. The ta- 
bles given in Carlson (1988, 1989, 1991, 1992) greatly 
simplify the reductions of the equations of motion to el- 
liptic integrals. The primary aim of this paper is to use 
Carlson's integrals to calculate all geodesic coordinates 
semi-analytically for the first time. 

Section 2 gives the geodesic equations in Kerr space- 
time. Sections 3 and 4 present the reductions to elliptic 
integrals and the specifics of our implementation. Sec- 
tion 5 outlines a variety of checks performed to ensure its 
validity and accuracy, and discusses the speed improve- 
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mcnt that should be expected from using an analytic 
code. Section 6 provides an overview of our code for 
readers who are not interested in all of its detail, and the 
code is applied to toy problems and test cases in Sec- 
tion 7. Finally, Section 8 discusses future work both in 
extending the code and in applying it to more realistic 
astrophysical situations. 



The signs of the integrals in r and 9 are independent 
and arbitrary, but are fixed for a given geodesic. It may 
seem odd that these equations lend themselves to the 
choice of r or ^ as independent variable to determine the 
cyclic coordinates t and ^. However, this is the natural 
outcome of the separation of the Hamilton-Jacobi equa- 
tion. 



2. GEODESIC EQUATIONS OF MOTION 

In Boyer-Lindquist coordinates {t,r,9,<p), the Kerr line 
element can be written. 



with the definitions 



-i2 /'„2 , „2n2 „2 A „;„2 



2 2,2 2 , 

p — r + a cos ( 



i:' = {r' + a')' -a'Asin'e, 



(1) 



(2) 
(3) 



where a is the angular momentum of the black hole and 
we use units with G = c = M = 1. 

Carter (1968) demonstrated the separability of the 
Hamilton-Jacobi equation for geodesies. 



^dS^ dS_dS_ 



(4) 



where S is Hamilton's principal function (the classical 
action) and A is an aSine parameter. The separation 
reduces the equations of motion to quadratures (Chan- 
drasekhar 1983) relating the coordinates r and 6: 



J y/R~ J 



7i' 



where 



R= [{r^ +a'^)E - aL^f 



e 



A [Q + {L, - aEf + 5ir2] 
Q- [ a2(5i -E'^) + Ll csc^ 6] cos" 6; 



(5) 



(6) 



(7) 



and the constants of the motion are the angular momen- 
tum about the black hole spin axis, L^, the energy, E, 
and Carter's constant Q. = (1) for null (timelike) 
geodesies. 

The equations of motion for the cyclic coordinates are 



t = \E + 2 j r[r^E - a{L^- aE)] 

r dr 

cf> = aj [{r'+a')E-aL,]^ + 



Ay/R 



/ 



{L^ csc^ - aE) 



7i' 



with 



cos^e 



d0. 



(8) 

(9) 
(10) 



3. REDUCTION TO CARLSON INTEGRALS 

In reducing the equations of motion from the previ- 
ous section, we follow closely the treatment given in Ap- 
pendix A of RB94. First change variables to (t, u, n, 

(j)) with p, = cos 9, u — 1/r. This set is more useful 
computationally, since the location of an observer at in- 
finity is mapped to u = 0. The domain of u is then 
< u < M+ < 1, where u+ is the location of the event 
horizon. Similarly, — 1 < /U < 1. Then the definitions 
= Q/E'^, I = Lz/E, and 7 = E/m put the equations 
of motion in dimensionless form. The integral equation 
relating u and is 



f f 



du 



where 

M 
U 



= + id^ -q' 
= (1 - 7-2) -h 27-2^-^(52 
2[{a-lf + q^]'i 
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(11) 



(12) 



(13) 



and = (1 — 'y~'^)a^. This paper only considers null 
geodesies, so that 7"^ = Q throughout. The arbitrary 
signs have been written explicitly, and are chosen to be 
Sx = sign(i:), where a dot refers to a derivative with 
respect to affine parameter. This is done so that both 
sides of (11) are always positive. The equations for the 
other coordinates become 



t-to 



+ 



f 2 2 diJ. 

2a{a - l)u^ + a^v? 



■"/ 



1 du 



v?{ulu+ - l){u/u- -l)s/U 



2(a-l)u + l 



du 



{u/u+ - 1){U/U- - 1) ^/u' 



(14) 



(15) 



where u± = [1 ± \/T~~a^]^^. The limits of integration 
have been omitted due to complications in accounting for 
turning points. This is discussed in more detail below. 

Given initial and final values of u and fi, we can com- 
pute t and (j). Since the /x integral is easier to invert 
and this method is of more general utility, u is taken 
as the independent variable and the goal is to solve for 
fif given /if). Mo and Uf. In certain applications it is 
more convenient to choose /U as the independent vari- 
able. For example, in the case of thin disk accretion we 
know the inclination angle as well as the value of fi where 
the geodesic intersects the disk. Section 4 gives solutions 
for Uf given /Uo; ^^'^ ^0 to handle these cases. 
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TABLE 1 
Reduction of I, 









l-^ T CI TVl O't" Oy" \-< CtTi <Tt^ 






RB94 


1 


Cubic (3 real) 


a = 0, 


+ > 27 or u< U2 


(-«1,1), (U2,-1), (W3,-1) 




1,3,8,10 






a ^ 0, g2 


= 0,\l\7^\a, 








2 


Cubic (3 real) 


a = 0, 


+ > 27 or M > n3 


(-Ml,l), (-M2,l), (-«3,1) 


2,4,9,11 






a ^ 0, g2 


= 0,|/|7^|a, 


{-ui,iy,{2ui[{a - If + q^]-\f /ui,l) 






3 


Cubic (1 real) 


a = 0,q2 


+ < 27 or a ^ 0, q2 = 0, Z ^ a 


[0,u+) 


5,7,12 


4 


No roots 


q^=0,l-- 


= a 




[0,«+) 


6 


5 


Quartic (2 real) 




> 


l),gs(«4, -l)'*;([-a2g2Mi«4]-i, [n-i + uj^]/, 1) 


[0,u+) 


13,19 


6 


Quartic (0 real) 


a ^ 0, ij2 


< 




[0,u+) 


14 


7 


Quartic (4 real) 


a 7^ 0, g2 


^ 0, tt < 1(2 


(-Ml,l), (m2,-1), (w3, -1),(«4, -1) 


[0,«2]'' ^ 


15,17 


8 


Quartic (4 real) 


a ^ 0, g2 


^ 0, U > U3 


(-Wl,l), (-W2,l), (-U3,1),(U4,-1) 


[U3,u+f 


16,18 



^hi is found from solving Eq. (21) and selecting one of the two real roots, d, e are defined in Eq. (22). 
^When U2 = the domain of u is [0,W2). 
^When U2 — lis, the domain of u is (u3,n+). 
"^98 = •siS"(9^). 



3.1. Reduction of lu 

Call the left-hand side (LHS) and right-hand side 
(RHS) of (11) Ifj, and /„ respectively, and start with the 
reduction of /„: 

Except in the special case with a = l,q^ = 0, U{u) is 
either a quartic or cubic and its roots arc denoted Ui with 
i = 1 — 3, 4, and ordered increasingly. If real, ui < and 
in the quartic case, U4 > 1 or U4 < 0. They arc of no 
physical significance. When all roots are real, the allowed 
regions for the integrand are u > U3 and w < M2 so that U 
is positive. Thus the roots are the turning points for null 
geodesies starting outside U2 and inside M3 respectively 
in both the cubic and quartic cases. There can be no 
more than one turning point, since the allowed region 
is bounded on one side either by infinity or the event 
horizon. When one or both pairs of roots are complex, 
there is no turning point in u. 

Upon encountering a turning point, the sign of u is 
reversed, so that the total integral is the sum of the inte- 
gral from uo to the turning point and that from u/ to the 
turning point. The idea is to ensure that the integrals 
in u and n monotonically increase along a geodesic. In 
a sense this allows the independent variable to take the 
place of the affine parameter, which cannot be used since 
it is a function of u and jj,. 

Carlson (1988, 1989) contain formulas for evaluation 
of integrals of the form 

[P]^ rUi^i + htr'/^dt; (17) 

-'y i=i 

with all quantities real, x > y, and cii + bit > for 
y < t < X. The form of a given integral is described 
by the vector [p], which contains the powers, Pi, of the 
factored roots. Cases with one or two pairs of complex 
roots are handled in Carlson (1991, 1992), where they 
are written in terms of real quantities as 



[p]= / {f+gt+ht^r/^ n {cii+hitr/Ht (18) 

•'y i=i,4,5 

for one pair of complex roots or 

2 

[p] - / n(,A + 9it + hity^'\a^ + hty-^'^dt (19) 
■'y 1=1 

for two. In using this form, it is assumed that each power 
Pi of an irreducible quadratic is written twice in the vec- 
tor [p]. In other words, when one pair of roots is com- 
plex, p2 = P3- When all roots are complex, p2 = P3 and 
Pi = Pi- 

To ensure that x > y in cases where a turning point 
may be present, integrals are written in pieces involving 

the relevant turning point, u,, and the number of turning 
points along the portion of the geodesic being followed, 
iV„ (either or 1): 




The Carlson papers reduce all elliptic forms to a set 

of four fundamental integrals, known as the R-functions 
(Press et al. 1992), which replace Legendre's integrals of 
the first, second and third kind. They arc all integrals 
from to 00 and hence don't require a limit of integration 
to be a turning point, greatly simplifying complex root 
cases where no physical turning point is present. This is 
one of many advantages of Carlson's approach. As is the 
case for Legendre's formulation, any elliptic integral can 
be reduced to a sum of Carlson's R-functions. Where 
Legendre integrals are used in this paper, they are cal- 
culated in terms of the R-functions using the formulas 
in Press et al. (1992). The integrals encountered in this 
paper are always of the form p = [— 1,— 1,— 1, —1,^5] for 
quartic cases and p = [— 1, — 1, — 1,^5] for cubic cases. 
Thus the form of coordinate integrals in the following 
will be specified by p^ alone. 

To maintain as much generality as possible, all inte- 
grals are written as above in terms of their roots. In cubic 
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cases the roots are found from solving the cubic equation, 
while for quartic cases they arc found numerically using 
the routine zroots.f from Press et al. (1992). Finally, 
instead of writing out the explicit formulas from Carl- 
son's papers and going through the algebra separately in 
each case, we have written routines for each case. This is 
much simpler and of more general utility, since numerous 
integrals must be done to calculate the coordinates of a 
point along a geodesic. 

The integral /„ has Ps = and is given by Carlson 
(1989) Eq. (2.12) for real roots for cubic cases. Quar- 
tic cases arc found in Carlson (1988) Eq. (2.13) for real 
roots and Carlson (1992) Eq. (2.36) for all complex roots. 
The quartic and cubic cases with a single pair of complex 
roots are given by Carlson (1989) Eq. (3.8). The nec- 
essary arguments to the Carlson routines are listed by 
case in Table 1, along with case definitions, appropriate 
domains of u, and the corresponding cases in Appendix 
A of RB94. 

As can be seen from Table 1, writing formulas in terms 
of the roots of U has the advantage of unifying many 
disparate cases from previous work. Equal roots cases, 
which describe orbits approaching the unstable circular 
photon orbits, cannot strictly speaking be treated iden- 
tically to other real roots cases as shown in the table. 
Here, integration to the turning point diverges. The code 
flags for these cases and integrates them directly from 
Uq to Uf, and the arguments listed in the table are still 
valid. In practice, however, except for the well known 
Schwarzschild unstable circular orbits with q"^ + l"^ = 27, 
equal roots cases are almost impossible to trigger. This 
is because the Carlson routines as written maintain ac- 
curacy until \u2 — < 10~^^, which is usually more 
precise than the determination of the imaginary parts of 
the roots. 

For one pair of complex roots, the arguments /, g 

and h are foimd by setting U{u) = qs{ui — u){u — 
ui){f + gu + hv?), where qs — s\gn{q^), and match- 
ing powers of u. When all roots are complex, setting 
U{u) = {fi+giu + hiu'^){f2+g2U+h2u'^) yields five non- 
linear equations for our six unknown coefficients. The 
degree of freedom is used to simplify the equations, and 
a sixth degree polynomial is solved numerically for hi: 



c 



2'-- 

e 



hf-hl i=hi+l 



where 



= 0, 

(21) 



c = a^ -l-" -q\ d = 2[{a-iy + qX e = -a^q^. (22) 

The only pair of real solutions to this equation corre- 
spond to the values of hi, /i2- 

As a full example of one of these reductions, consider 
case 5 from Table 1 with uq < v f (.s„ = 1). This is 
the Kerr case with no physical turning points. From 18, 
we see that 61 = 1, 64 = —qs, ai = — wi, 04 = qsU4, 
X = Uj\ y = 2iQ. The sign qs is used to keep each factor 
positive. Matching the powers of U (u) as described above 
gives / = -qs/{uiUie), g = {m + ui) / {U1U4) f , h = 1. 
Following Carlson (1991), we define 



C = ^Jf^gx + hx-^, ij ^ Vf + gy + hy^, (24) 
Cij = \J 2fbibj — g{aibj + ajbi) + 2haiaj, (25) 

M= ^''^^^ + '^^''^V (€ + ^)^-M^-^)^ (26) 
x y 

LX=M'^ + ch±CllCiA. 



Then, 



Rf{M\lI,LI). 



(27) 



(28) 



Rp is computed using the routine from Press et al. 
(1992). Equations for Carlson elliptic integrals with 
P5 can similarly be found in the Carlson papers 
listed above. 

3.2. Inversion of I^j, 

Next, the /,j integral needs to be inverted to solve for 
lif. As with U{u), the roots of the biquadratic M{ii), 
M±, determine the physical turning points in ^. When 
M_ > 0, there are four real roots and the orbit cannot 
cross the equatorial plane. The physical turning points 
correspond to the two roots with the same sign as /zq 
and are denoted /u± = sign.{iJ,o)y/M±. When M_ < 0, 



the physical turning points are /U± = ±1 



and are 



symmetric about the equatorial plane. We can calculate 

the number of times the geodesic has crossed a /j, turning 
point from the magnitude of the /„ integral. This is done 
by noting that the maximum value of fj^^ is J^"*" and its 

minimum value is zero. In this derivation the integrand 
d^/^fM, common to all integrals, is omitted. Then, for 
Su. = 1, 



/ +{N-1) <Iu< +N , 



(29) 



where N is the number of turning points reached in /j,, 
and fj,± are the upper and lower turning points in ^. 
The integrals are written in these pieces so that they 

arc always positive, as required for use with Carlson's 
integrals. This condition can be written more concisely 
as 



N 



(30) 



where [] is the ceiling function. If = —1, then the 
first turning point reached is /i_ . The condition can then 
be written 



+ (^-1) 



Mo 









r 


<Iu<- 











y, 



(23) 



(31) 

Using f'^ = f^^ — f"*", we can rewrite this in terms of 
the same integrals used above: 

- +N <Iu<- +(A^+1)/ • (32) 
Finally, 



Analjrtic Kerr Photon Orbits 



5 



N = 



I + 



(33) 



and [ J is the floor function. To write out the general so- 
lution for /„ = I/j, for arbitrary number of turning points 
and Sfj,, we include coefficients for the various pieces of 
the I/j, integral: 



rf^+ ri^f ri^+ 

Iu=ai / +a2 / +a3 / 

J Lin J U— J U- 



(34) 



The coefficients arc functions of and N determined 
by writing down specific cases. For example, tti reflects 
whether the integration is positive or negative from hq 
to and is easily seen to be ai = s^. Similarly, a2 
reflects whether the last turning point reached is /i_ or 
Thus the coefficient is a2 = s^(— 1)^. The third 
coefficient is slightly more complicated and turns out to 
be 



as = 2 



2N + 3-S. 



- 1. 



(35) 



Armed with the number of turning points and the co- 
efficients, wc solve for /U/ by inverting the second integral 
on the RHS of (34): 



J U.— 



M Q!2 



1 / ft^+ rf^+ 

lu- cti / -as 



(36) 



Calling the RHS / and writing out the square root on 
the LHS for the general case (o j^O, 0) gives 



1 



\a\J^_ V(M+^7?)aI^^Ml)' 



(37) 



Carlson (2005) contains a table for inverting integrals 
of the form 



dt 



(38) 



where all quantities are real, x>y,0<y<x and either 
y = 0, a; = oo or one limit is a root of the integrand. The 
latter case applies here. 

3.2.1. M_ > 

When M_ > 0, all requirements are met as written, 
and 



jif = H-nd{J,k), J = /i+|a|/, fc^ = 1 — (39) 

where nd (J, k) = l/dn{J, k) and dn is a Jacobi-EUiptic 
function. The fi integral terms in / are calculated as 



A*/ 



dji 



1 



, = -rF{x,k) 



(40) 



where F{x, k) is Legendre's integral of the first kind 
(Abramowitz & Stegun 1965), x = j^JTjvf" -, A = |a|/i+ 

and k is the same as above. The integral between turning 
points is the complete elliptic integral K{k). 



3.2.2. M_ < 

When M_ < 0, y < in (38) so that (39) is no longer 
valid. Since the integrand is an even function of we 
can write 



1. /''^ 



(41) 



which is in the correct form, except that —fif can be 
negative. This causes no problems. In this case. 



Ml 



Hf = H-cn{J,k), J = ^ M_\a\I, k = 2 _ ' 

(42)" 

and we've used /i_ = — /x+ for M_ < 0. The fj, terms in I 
are computed the same as in (40), with k defined in (39), 

X = ^1 - ^ and A = \a\-s/M+ - M_. The integral 

between the turning points here is twice the complete 
elliptic integral K{k). 

3.2.3. g2 = 

A special case is encountered when = 0. M(/Lt) has 

a double root at /i = 0, causing /^^ to diverge there, 
and preventing these orbits from reaching the equatorial 
plane. Hence, they have at most one physical turning 
point. In this case is elementary, and the solution for 
M/ is 

M/ = M+ sech [|aA<+|/„ - s^si sech"^(/xo/M+)] , (43) 
where Si = sign(/io)- 

3.2.4. a = 

Finally, when a = (the Schwarzschild case) the fif 
integral is again elementary. The solution for /x/ is then. 



jjLf — H- cos 




ai cos 



0:3 TT 



(44) 



3.3. t and 



coordinate integrals 

Given the solution for /x/, equations for the coordinates 

t and can be reduced to elliptic integrals as well. Each 
coordinate is expressed as a sum of integrals over u and 
/i. As is done above, the u terms are reduced to Carlson's 
formulation, and the /x terms to Legendre's. 

The n integral term in (14), which we'll denote T^, can 
be written as a single Legendre integral of the 2nd kind. 
For example, the /io term in the M_ < case is reduced 
as follows: 



7; = |a| 



^(mI-/x2)(/,2_M_) 



= \a\n+ 



dt 



^(l-t^){l-t^-^) 



l-t' 



dt- 



1^ 



Io ^(l-f2)(l_i2_ 

= AE{x,k)+a'^M_Iy,, 



M_ 



) 



+ OfM_Iy, 



(45) 
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where E{x, k) is the Lcgendrc integral of the second kind 
with arguments x and k defined in the previous section. 
The substitution t = a/I — no/ ii+ is made between hnes 
one and two, and M_//x^ is added and subtracted from 
the numerator between lines two and three. In the > 
case, is given by the first term of the above formula, 
with the arguments A, k, x for that case given with the 
solution for /i/ in Subsection 3.2.1. 

The fji term in the (p component formula (15) can be re- 
duced to a Legendre integral of the 3rd kind in analogous 
fashion. For the M_ < case we proceed as follows: 



-llu 



-llu + 



\a\iJ.+ 



dt 



The first term is from T„ and the second term is T^. 

Component integrals are calculated the same way as 
J„ or 1^ respectively. That is, component integrals are 
calculated in pieces using the appropriate coefficients as 
described above while u component integrals are calcu- 
lated with reference to the physical turning point, if one 
exists. These are all the integrals required to compute 
null geodesies in Kerr spacetime. These equations for the 
0, t coordinates are written in Boyer-Lindquist coordi- 
nates. For certain applications, Kerr-Schild coordinates 
are used instead. We note here for completeness the ana- 
lytic transformations between our Boyer-Lindquist coor- 
dinates and these Kerr-Schild coordinates {i,u,jl,<p) (Font 
et al. 1999), 



i = i + logA-w^log| ^ ^ ^7^= K (50) 



^+ + M+^i^ ^(l-t2)((l_^)_t2) 



l-u[l- VI - a2] 



-llu + 



I 



-n(n; X, k) 



A{1 - M+) 

where n(n; x, k) is the Legendre integral of the 3rd kind 

2 

and n = ^ ■ The formula for the M_ > case is the 



(46) 



1 = 4> — aur log 



1 - it[l + Vl - a2 



l-u[l- VI - a?] 



same, with n 



u = u, 



(51) 

(52) 
(53) 



The formula for the M. 

1-M+ 



and the other arguments defined 
in Subsection 3.2.1 above. Note that we are using the 
sign convention for n from Press et al. (1992), which is 
opposite that in Abramowitz & Stegun (1965). 

T„, the u integral term in (14), is expanded with partial 
fractions, and after a little algebra is written 



Tu 



2a(a -l)-{ h 

u+ 



1 



du 



(2a(a-0 + ^ + ^)/ 



{u/u+ -1)VU 
1 du 



(u/u- - 1) Vf7 



+ 



it U\) J U\/U Ur J 



du 



uWU 



(47) 



where Ur 



-(2vr 



is negative. Three 



U+U- 

U^—U — 

of the terms have — —2 and one has ~ —4. When 
a limit of integration is at infinity {u = 0), this integral 
blows up, as it should. In practice, the code picks a non- 
infinite starting radius large enough that the geodesic 
trajectories from infinity to the starting radius differ neg- 
ligibly from their flat space counterparts. 
Then, 



SiiUf 



— + 2(a- 
u+ 



I) 



1 



du 



±+2{a-l)] f / 

u- J J (u/u^ - 1 



{u/u+ -1)VU 
du 
If 



,(48) 



where both integrals are already calculated as part of T„. 

Finally, the dimensionless affine parameter can also be 
calculated along the path from (10) without any addi- 
tional integrals: 



The transformations are valid outside the event horizon, 
where A and the numerator of the other log terms are 
positive. 

4. SOLUTION FOR Up 

For some applications, it is preferable to use /it as the 
independent variable and solve for Uf given uq. In par- 
ticular, consider geodesies connecting an observer at in- 
finity with a thin, equatorial accretion disk. The initial 
polar angle is the inclination of the observer. The final 
polar angle is 7r/2 (/x/ = 0), and we solve for the ra- 
dial coordinate where the ray intersects the disk. This 
method, however, is of less general utility than that de- 
scribed above. Even in simple geometries, the number 
of turning points in fj, along a geodesic is not known in 
advance as it must be to use n as the independent vari- 
able. One way around this is to calculate all geodesies 
connecting the observer with the disk for a fixed num- 
ber of /i turning points (Cunningham & Bardeen 1973; 
Viergutz 1993). 

The approach in solving for w/ is the same as in solv- 
ing for Uf. The integral is computed as a Legendre 
integral of the first kind. Given the number of turning 
points, is computed in pieces as shown above using 
the coefficients 0:1,2,3- 

After finding 7^, we invert /„. This inversion 
ranges from relatively straightforward to algebraically 
formidable. As examples, we discuss cubic and quar- 
tic real roots cases in detail. Table 2 gives the solution 
for Uf in all cases. This problem was first addressed by 
Agol (1997) and the solutions here are from its Table 5.2 
with some modification. 

For our first example, consider the first two cases of 
Table 1 where there are three real roots. The integral to 
invert is 



f du 2 f 



(x^dfj. 



(49) 



\Jua 



du 



± 



du 



Uf 



(54) 
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TABLE 2 
Solution for uf 



a 


f 












m\ 


Cl 


C2 


C3 



''I + U2icd?J 
ui + u:iidc^J 

C2+ui — (c2—ui)cnJ 
l-\-cnJ 

•I.- " 1- ('J,- " I'Ti - "I'-.l 



(C4 -g,,c-, )C71 J + q,, C4+C5 
n(l + c|)^e.J 



1 — C2 sc J 

1 — C2Sn^ J 
1 — C2S7l-^ J 



Cl[-fM - -f«('"0,«2)] 
Cl[V + Iuiu3,Uo)] 
Cl[I^ + Iu(ui,Uo)] 

ci[/m + Iu{ut,uo)] 

SuCl [Ifi + Iu{cs,Uo)] 

cill^i - /«(no,«2)] 
Cl[I^ + Iu(u3,Uo)] 



"31 

"31 
1_ . 6"i+C3 

O 



8C2 



(c4+<;.s'-.) 



4C4C5 

/ e4-e5 \ ^ 
VC4+C5/ 
"41 "32 
"42"31 
"41 "32 
"42"31 



^(C4 + C5) 

2 
2 



V«l(3«l +C3) 2±| 



— (C4 —05)-^ 
(c4+c5)^-4n^ 

"21 
"31 
"43 
"42 



M42U31 
«42«31 



^snj — S7i{j. mi), cnJ — cn{J^ mi), scJ — sn{J^ nii) / cn{J, mi) and mi — 1 — A: is used instead of A:. 
^qs — sign(g^). If q.y — 1 then Ua = U4, Ub = wi. Otherwise, Ua = wi,Ub = U4. 



*^ Complex roots are written as m ± in, p ± ir and are ordered so that m > p and n > 0. 



TABLE 3 

AuxiLLARY Constants used in table 2 



£4 C6 

6 V^CyTr^-fj)^~+~(nr+~r)^ ^ (m — p)^ + (n — r)^ 



where u+ is the relevant turning point: U2 or U3. Denote 
the first term , and write the second term in terms of 
the roots of the integrand: 



± 



Vd Jut 



du 



'us \/{u- Ui){u - U2){U - U3) ' 

(55) 

where d = 2[{a — /)^ + q^]. This can be put in the form 
(38) with the substitution z = — u\: 



J = V^I^^/, fc2 = ^^H_Jfl, (60) 



dc(J,fc) = ^4TTT' cd(J,fc) 



U3 - Ul 

_ cn(J, k) 



cn{J, k) dn{J, k) 

where cn, dn are Jacobi-Elliptic functions. Note that the 
result does not depend on whether or not a turning point 
has been reached, since both cn and dn are even in J. 
When U{u) has four real roots and u <U2, 



/■V'M+-«l 
J ,/uf-Ul 



±1 



dz 



^Uf-ui \/{z'^ + (Ul - 'W2)][-2^ + («! - Ws)] ' 



where 



(56) 



(57) 



and lu^ is determined from the same Carlson formulas 
as for /„ above. 

Comparing (56) with (38), we see that a\ = u\ — U2, 
a2 = u\ — U3 and 61 = 62 = 1. If m+ = wa, then the 
limits of integration must be switched, since by definition 
X > y. These integrals correspond to the third row, third 
and fourth columns of Table 1 from Carlson (2005), and 
the solutions for Uf are 



with 



Uf = ui + {u2 — Ul) cd'^{J, k), 
= ui + (us - Ul) dc^(J, k), 



Uq < U2 (58) 
uq > U3, (59) 



I -I =+— [ 
" Vein 



du 



Uf \/{u- Ui){u - U2){u - U3)(W4 - u) ' 

(61) 



where e = \aq\. With the substitution z = y this 
becomes 



where 



^JXz"^ + (Wi - U2)\\z'^ + (Ul - Us)] ' 

(62) 



(63) 



Again comparing with (38) and using Carlson (2005), 
we find 

Us (U2 - Ul) Sn^ - U2(WS - Wl) 

"/ = —r:. -. (64) 



(u2 - ui)sn2 - (us - Ul) 



where 



8 



Dexter and Agol 



sn = sn{J,k), J = \J {u\ — U'i){ui — u\)l (65) 

(■U4 - M2)(M3 - Wl) 

Again the result is independent of whether or not a 
turning point is present. In (54) above, the sign of the 
second term on the RHS depends on whether a turning 
point is present. This allows us to determine the number 
of turning points in u. 

When complex roots are present, the reduction to 
standard form (38) is much more difficult. It is dis- 
cussed in Erdelyi et al. (1981), and relevant formulas for 
the inversion can be found there and in Byrd & Fried- 
man (1971). In particular, our cases 3,5 are from Byrd 
& Friedman (1971) equations 239.00 (p86) and 259.00, 
260.00 (pl33,135). Our formula for case 6 is based on 
Erdelyi et al. (1981) Table 2, p310-311. The intricacy of 
these reductions demonstrates the advantage of Carlson's 
method. The computation of integrals is equally efficient 
with complex or real roots. Unfortunately, when inver- 
sion is required, Carlson (2005) is only a somewhat more 
compact version of Legendre's original notation and of- 
fers no real advantage over previous work. 

5. CODE CHECKS AND SPEED TESTS 

Using the solution for /iy, the equality of and 
has been checked to machine accuracy (at least 14 sig- 
nificant digits in all cases). Once found, can be used 
as an input to recover u/. In this way, the two routines 
have been shown to agree in all cases. Precision in the 
calciilation is limited by error in the determination of 
the roots of U{u). An advantage of using Carlson's for- 
mulation is that all component integrals are computed 
without reference to the complex roots, which often have 
less numerical precision than real ones. Formulas for u/, 
discussed in Section 4, are also written in terms of real 
quantities leading to higher accuracy. 

Certain special cases can be integrated analytically for 
all components, providing independent checks on com- 
ponent integral formulas. These include /U cases with 
= 0, u cases with = 0, I = a and u cases with 
equal physical real roots, corresponding to unstable cir- 
cular photon orbits. In all of these cases, the component 
integral formulas above agree with the analytic results 
to machine accuracy. The component integral formu- 
las above also reduce to those derived separately for the 
Schwarzschild case, a = 0. All of the formulas given 
here for /x/ agree with those in tables 1-2 of RB94. The 
Schwarzschild formulas were also tested against the ap- 
proximate formulas given by Beloborodov (2002). 

Further, the implementations of Carlson's integral ta- 
bles have been checked extensively using the Mathemat- 
ica MIntegrate function. The same is true of the t and 
(j) formulas, as well as the individual integral components 
7„, Tu, and 7^, T^, 

The R-function routines maintain accuracy until a < 
10~^ or < 10^^°. If such parameters are encountered, 
the code will give a warning and set the offending value 
to zero. 

The geodesic computations have been checked against 
calculations done by the code used in Falcke et al. (2000) 
and are found to be in excellent agreement. The FOR- 
TRAN implementation of our code is found to be faster 



than that one by a factor of about 5, due to the fact that 
our code computes the minimum number of R-functions 
possible, and shares them between routines when neces- 
sary. The code from Agol (1997) was found to be ~ 100 
times faster than numerical integration in the case of 
tracing geodesies from infinity to a thin disk. This is an 
optimal problem for an analytic code, since we can solve 
for the point where /i = 0, whereas a numerical code 
must integrate the geodesic from infinity until it reaches 
that point, and then zoom in on the intersection to find 
an approximate solution to the desired accuracy. In addi- 
tion, this example didn't include the t and (j) coordinates, 
which are sped up by a much smaller factor than u and 

As a lower bound for the speed improvement of our 
code over numerical integration, a routine was written 
to integrate the photon four momentum for all coor- 
dinates with respect to affine parameter using the im- 
plementation of the Bulirsch-Stoer method from Press 
et al. (1992). We then compared the integration of many 
points along a single geodesic starting from infinity with 
this numerical code and our analytic one. This is the 
ideal case for numerical integration, since the intermedi- 
ate points calculated along the ray are no longer wasted 
as in the first example. For the case considered with no 
turning points in u or /z, the analytic code was found to 
be faster by a factor of ^ 3. 

However, our numerical code for integrating geodesies 
is much simpler than a complete code would have to be. 
It cannot handle turning points, and requires knowledge 
of the afiine parameter on the ray in order to know where 
the region of interest in the integration is. In practice, 
turning points would have to be detected and a scheme 
for determining the region of interest in affine parameter 
implemented. Alternatively, a somewhat more compli- 
cated scheme such as the Hamiltonian method described 
in Schnittman & Bertschingcr (2004a) could be adopted. 
In any case, these additions would slow down geodesic 
computation. We conservatively estimate, then, that the 
lower bound for the speed advantage of our analytic code 
over numerical integration is a factor of ~ 5. A similar 
test for the case of integrating down to the thin disk 
yielded a speed difference of a factor of ^ 300, leading to 
an upper bound on the speed advantage of ^ 500, which 
is in good agreement with the naive estimate of multi- 
plying the speedup found by Agol (1997) by the speedup 
factor between our code and the one presented there. 
Then the range of speedup that can be expected by us- 
ing the analytic code described here is a factor between 
5 — 500 depending mostly on the application, but also on 
the specific implementation of the numerical integration 
code. 

In addition to being faster, the analytic formulation is 
much more flexible. It can calculate an arbitrary num- 
ber of points beginning and ending anywhere on any 
geodesic, provided that the constants of the motion can 
be calculated. This is exploited in the thin disk toy mod- 
els below, where we solve for the point Hf = 0. It could 
also allow, for example, a calculation of Compton scat- 
tering by tracing rays out from every point on a geodesic, 
and computing the scattered intensity into that point as 
a separate ray tracing computation. In any event, the 
flexibility inherent to an analytic method could allow 
for more sophisticated calculations in the future, which 
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Fig. 1. — Change in time vs. radial coordinate in the 
Schwarzschild metric for geodesies near the circular photon orbit 
(dashed line), as described in Section 6. 
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Fig. 2. Image of a near extreme (a = .998) Kerr black hole 
viewed from the equatorial plane. Images intensities are taken to 
be the afline parameter evaluated upon termination at the black 
hole or after returning to the starting radius. Intensities are scaled 
linearly from the minimum value outside the shadow to the maxi- 
mum. 

wouldn't be possible with a numerical code. The main 
disadvantage of using an analytic code is that the affine 
parameter cannot be used as an independent variable, 
which may be desirable for adaptive integration tech- 
niques in radiative transfer applications, for example. 

6. IMPLEMENTATION 

This section provides an overview of the various rou- 
tines used by the code described above, and examples of 
their use. The README file online covers everything in 
this section in greater detail. The FORTRAN 77 source 
file geokerr . f contains the main program as well as the 
key routines, geokerr, geomu, geor and geophitime, 
and supporting functions. Inputs are given through com- 
mand line prompt or a text file. Inputs from previ- 
ous command line runs may be saved for future use. 
These inputs include constants of motion for the desired 
geodesies, initial and final u and initial /i, the number 
of turning points in u (ignored if the constants do not 
admit physical turning points), and the sign of ii and fi. 



Constants of the motion are required and may be speci- 
fied either as the impact parameters at infinity (a, as 
is most convenient in ray tracing applications, or as the 
dimensionless angular momentum, /, and Carter's con- 
stant, g^. When any other information is not provided, 
the program assumes geodesies which trace out the entire 
domain of u, from the starting point and back or until 
the event horizon is reached. 

The program calls the main subroutine, geokerr, 
which calls geomu to fill in missing inputs and calcu- 
late /!/. Alternatively, geokerr can solve for Uf using 
geor. Subsequently, geophitime calculates the (f> and t 
integrals using the Carlson routines. The program loops 
over constants of the motion for a chosen initial polar an- 
gle and black hole spin. Results are written to standard 
output by default, and should be redirected from the ter- 
minal to a text file in most cases. It is also possible to 
input the name of the desired output file. The subroutine 
geokerr encapsulates most of the code functionality, and 
can fairly easily be adapted to another front end other 
than the program used here. See the README file ac- 
companying the code or online for more detail. 

As an example use of the code, consider tracing rays 
over a rectangular grid in —4 < a < 8, —6 < /? < 6 for 
a near extreme black hole, a. = .998. The observer is at 
infinity in the equatorial plane (/zq = 0), and 20 rays will 
be traced over each dimension. The input file for this 
situation can be found online.^ 

Output is arranged as follows. The constants of the 
motion are listed for each geodesic in the top line, fol- 
lowed by columns giving Uf, fif, At, Acj), A. The format 
used for output can be changed with a tiny modification 
to the source code. See the README file for more de- 
tails. Plotting the affine parameter evaluated at either 
the event horizon, or once the geodesic returns to its ini- 
tial radius, as a function of impact parameters for this 
data with 160, 000 geodesies produces Fig. 2 as explained 
below. 

For less standard batch runs, it may be necessary to 
generate the input file from a simple program. Consider 
a set of geodesies in the Schwarzschild metric (a = 0) 
to study the unstable circular photon orbits. The given 
parameters are chosen to be uq = 1/30, Uf — u+ = .5, 
/io = .9, (3 = and an array of values for a near \/27. 

The piece of code to write an appropriate input file 
is available online.^ Plotting the change in time as a 
function of final radial coordinate produces Fig. 1. 

7. APPLICATIONS/ VALIDATION 

We next describe a couple of relatively simple appli- 
cations of the code to ray tracing problems as further 
validation and as examples of its utility. The first is the 
simplest illustration of the black hole shadow, which tests 
the determination of the roots of U (u) and qualitatively 
parts of the time integral. Next are examples from the 
standard model of thin disk accretion. The disk image 
and simple spectrum from line emission test the routine 
that solves for Uf. The projection of a uniform grid at 
infinity onto the equatorial plane of the black hole also 
tests the calculation of (j), and hot spot emission pro- 
vides a time-dependent test. Finally, spectra and images 

^ http ; / / www. astro .Washington . edu/agol / geokerr / exfiles / abgrid . in 
^ http: / / www.astro.washington.edu/agol /geokerr / exfiles /inputex.f 
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Fig. 3. — Projection of a uniform Cartesian grid in the image 
plane to the equatorial plane of the black hole for /iq = I (top) 
and fiQ = .5 (bottom). Black hole spin is a = (left) and a = .95 
(right) , and the area inside the horizon is removed from each image. 
Compare to Fig. 2 of Schnittman & Bertschinger (2004b). 

of synchrotron radiation from spherical accretion quan- 
titatively test our radiative transfer routines. 

Ray tracing utilizes the simple relationship between 
points on an observer's instrument and the constants 
of motion of null geodesies. Consider the photographic 
plate at infinity as a function of the impact parameters 
a, (3, perpendicular and parallel to the black hole spin 
axis respectively. Images can be created by tracing rays 
backwards from points on the plate to the black hole. 
The parameters a, /? are easily expressed in terms oiq^,l 
using (Cunningham & Bardeen 1973) 

l = ~a{l-^,lf^ (67) 
q''=0'+nl{a^-~a''), (68) 

so that each point on the observer's photographic plate 
corresponds to a unique geodesic. 

7.1. Image in Affine Parameter 

As a first application of ray tracing, we can deter- 
mine the appearance of the simplest possible black hole 
shadow. The image "intensities" are taken to be the 
affine parameter evaluated at the termination of the 
geodesic-either when it terminates at the black hole or 
reaches a turning point and re-emerges to the starting 
radius. Affine parameter is a good proxy for the emis- 
sion in this case, since it is related to the proper length 
along a geodesic, which would be the observed inten- 
sity for constant emissivity and neglecting absorption. 
The dimensionless affine parameter. A', is given by (49). 
The equatorial plane result for a Kerr black hole with 
a = .998, to be compared to Bardeen (1973) Fig. 6, is 
shown in Fig. 2. The image shown here is 400 x 400. 

7.2. Thin Disk Accretion 

The next set of applications imagine the emitting 
source as an infinitesimally thin disk in the equato- 
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Fig. 4. — Image of an optically thick standard relativistic accre- 
tion disk around a near extremal black hole (a=.998). The disk 
has outer radius rout = 18, and the observer's inclination is 85°. 
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Fig. 5. — Normalized spectra of line emission from a thin ac- 
cretion disk at an inclination of 30° for various black hole spins. 
The emissivity is taken to be proportional to between the 

marginally stable orbit and Rout = 15. Compare to Schnittman & 
Bertschinger (2004b) Fig. 3. 

rial plane of the black hole (e.g.. Page & Thorne 1974, 
Shakura & Sunyaev 1973). 

7.2.1. Grid Projection 

The first check of the code for this case is in visualiz- 
ing the projection of a uniform grid at infinity onto the 
equatorial plane of the black hole. This is done by solving 
for the final radius, u/, and azimuth where the geodesic 
intersects ^/ = 0. Then, the new grid points are calcu- 
lated using pseudo-Cartesian coordinates (Schnittman & 
Bertschinger 2004b): 

x= ^/r"^ + 0? cos (j), y = \/ r"^ + sin (p. (69) 

The result of this projection for two different initial 
observer inclinations and black hole spins is shown in Fig. 
3, and agrees with Fig. 2 of Schnittman & Bertschinger 
(2004b). The gravitational lensing effect can be seen in 
the pictures with /xq = .5 as the bunching of grid points 
behind the black hole, while frame dragging is evident in 
those with a = .95 

7.2.2. Thermal Disk Images 

As a next step, we can use the standard thin disk re- 
sults for the radial temperature profile (e.g., Krolik 1998) 
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to produce images of the disk at various inclinations as- 
suming it is optically thick everywhere, so that the inten- 
sity is that of a blackbody. Finding the radii of emission 
from a grid in impact parameters and calculating the in- 
tensity at each of these points produces an image of the 
disk as seen by a distant observer. The result for an in- 
clination of 85° and black hole spin a — .998 is shown in 
Fig. 4. The image shows the effects of relativistic beam- 
ing of the emission from gas moving towards the observer 
versus the redshift of that moving away, as well as the 
bending of the light from gas behind the black hole. 

7.2.3. Line Emission 

Next, following Schnittman & Bertschinger (2004b) 
and Bromley et al. (1997) we consider monochromatic 
emission from the disk, and give it an inner (outer) ra- 
dius, Rin = Rms {Rout = 15), where Rms is the loca- 
tion of the marginally stable circular orbit (e.g.. Page & 
Thorne 1974). The emissivity is weighted by uj, physi- 
cally motivated by the fact that we expect the tempera- 
ture of gas in the disk to increase with decreasing radius. 
The observed intensity is computed by exploiting the in- 
variance of Ii./v^ (Misner et al. 1973), 



(70) 

where g = vq/v is the redshift, and (u) is the ob- 
served (emitted) frequency. To see the effect of black 
hole spin on the emission in this c;ase, we calculate 1,^^ as 
a function of g for several values of a by calculating the 
intensity of rays at a location with redshift in a certain 
range of g, and integrating them over the photographic 
plate. The result is plotted in Fig. 5, and is in excel- 
lent agreement with Fig. 3 of Schnittman & Bertschinger 
(2004b). At higher black hole spin, the marginally stable 
orbit is much closer to the black hole where the redshift 
is much stronger, leading to a higher relative magnitude 
and broadening of the low frequency peak ("red wing"). 

7.2.4. Rotating Hot Spot 

Finally, to test the time-dependence of the code, con- 
sider a circular hot spot of finite radius Rspot = -5 orbit- 
ing in the equatorial plane of a Schwarzschild black hole 
at its marginally stable radius {Rms = 6). The emis- 
sivity of the spot is taken to be Gaussian in the locally 
flat space near the hot spot (for details, see Schnittman 
2006), 



j{x) (X exp 



^spot 



{t)\' 



9R2 



(71) 



where j is the monochromatic emissivity. For some ob- 
server coordinate time, t, the time delay and azimuthal 
position from the observer to points on the disk are used 
to determine where on the photographic plate the sepa- 
ration between geodesic and hotspot are less than ARspot- 
For these points, the Gaussian emissivity and observed 
frequency (redshift) are tabulated. Repeating this proce- 
dure over a period of the motion gives a time-dependent 
spectrum, which is shown in Fig. 6 for an observer incli- 
nation of 60° (/^o = -5). This figure is in good agreement 
with Fig. 4 of Schnittman & Bertschinger (2004b). 

Integrating over frequency (redshift), or equivalently 
over the impact parameters, gives the light curve. Fig. 7 
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Fig. 6. — Spectrogram of a circular hot spot of radius Rspot = -5 
at the marginally stable orbit of a Schwarzschild bla<;k hole. The 
observer is inclined at 6o = 60° . Compare to Fig. 4 of Schnittman 
& Bertschinger (2004b). 
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Fig. 7. — Light curves of the hot spot described in Fig. 6 for var- 
ious inclination angles. Intensities are normalized individually to 
the integrated intensity over each orbit and scaled to the maximum 
intensity from all inclinations. Compare to Fig. 6 of Schnittman 
& Bertschinger (2004b). 

shows the light curves of the hotspot for several inclina- 
tion angles. As the observer approaches edge-on viewing, 
the light curve becomes sharply peaked by a combination 
of the Doppler beaming of the spot as it moves toward 
the observer and the large gravitational lensing of the 
spot as it goes behind the black hole. The plot here is 
in excellent agreement with Schnittman & Bertschinger 
(2004b). 

7.3. Radiative Transfer 

In more realistic astrophysical applications, the source 
is not a delta function at a given inclination, and the 
intensity along a ray can be written more generally as 
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Fig. 8. — The spectrum of synchrotron radiation from optically 
thin spherical accretion onto a stellar mass black hole. The solid 
line is the ray tracing result, and the plotted points are the analytic 
results. The two curves agree to within 5% at low frequencies, 
where the radiation originates at larger radii and the bending of 
light should be unimportant. 



J ray ^ 
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If absorption can be neglected, dl^, = ji,dl where 
dl = —pau"'dX is the proper length differential measured 
along the ray, is the photon four-momentum, is the 
four-velocity of the emitting particle and A is an afBne 
parameter. 

The observed intensity is then 



'-dX 



(73) 



where ji, is the emission coefficient in the rest frame of 
the gas, and A is now the dimensionless affine parameter 
used above. It is calculated from (49) and used as the 
independent variable along the ray. When absorption is 
included, the solution to the radiative transfer equation 
between affine parameters Aq and A reads (Fuerst & Wu 
2004) 

hoW = <7'V.(Ao)e---(^«) + t e-(-''(^')--''(^°))32j, d\' , 

(74) 

where r,y = J a^dl is the optical depth. Throughout this 
paper we neglect scattering contributions to the emission 
and absorption coefficients. 

7.4. Synchrotron Radiation from Spherical Accretion 

The code described above in conjunction with a routine 
to perform radiative transfer along rays is now applied to 
the particularly simple case of a stellar mass black hole 
at rest with respect to the interstellar medium with a 
temperature at infinity of IC* K and a density at infinity 
oflcm^''. Ionized hydrog 



i.cc;rctes onto the black hole, 
and the magnetic field threading the gas effectively cre- 
ates collisions, so that the accreting gas can be considered 
a perfect fluid. In the model, magnetic turbulence es- 
tablishes an equipartition of magnetic and gravitational 
energy (Zeldovich & Novikov 1971). Then 



Fig. 9. — Spectrum of synchrotron radiation from spherical accre- 
tion onto a stellar mass black hole. The solid line is the ray tracing 
result including absorption. The spectrum is heavily attenuated at 
z/() < 10"'^ Hz, and in this region follows the optically thick approx- 
imation of thermal emission from the r = 1 surface. The spectrum 
agrees well with the emission only model for i/q > 10^^ Hz 




E 0.4 



Fig. 10. — Image of a spherically accreting Schwarzschild black 
hole aX V = 10^^ Hz as a contour plot and a 1-d profile. 
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and cgs units are most convenient in the analytic calcu- 
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lation. We assume an adiabatic equation of state with a 
piecewise adiabatic index (Shapiro & Teukolsky 1983), 
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(76) 



where rup, rUe are the proton and electron mass and T 
is the temperature in units of proton rest energy. Then 
the fluid equations arc non-linear, and can be solved nu- 
merically (Michel 1972) to find the temperature and fluid 
velocity as functions of coordinate radius. 

The dominant form of radiation produced is syn- 
chrotron radiation from the inner part of the accreting 
sphere, where the electrons are ultrarclativistic (Shapiro 
1973b). In this case, the emissivity can be well approxi- 
mated analytically. Shapiro (1973a) performed the rela- 
tivistic radiative transfer by approximating the photons 
as traveling on null geodesies in Minkowski spacetime, 
and calculating gravitational redshifts as well as the pho- 
ton Doppler shifts along these paths. 

Shapiro's formula for the radiated spectrum is 
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where v{r) is the proper velocity seen by a stationary 
observer and 
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is the critical angle at which the light is recaptured by 
the black hole. 

The synchrotron emissivity for thermal, ultrarclativis- 
tic electrons averaged over polarization and solid angle 
assuming isotropic emission in the rest frame is given by 
(Pacholczyk 1970), 
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and where 



F{x) 



K5/3{y)dy 



(81) 



(82) 



is the synchrotron function. Mahadevan et al. (1996) 
have approximated I{x) above analytically by matching 
the asymptotic forms for large and small x. They find 
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(83) 

Note that this function is denoted I'{x) by Mahadevan 
et al. (1996), and has a maximum error of « 2.7%. The 
spectrum is calculated by integrating (77) numerically. 

To compare with these results, the ray tracing code 
is used to create an image of the synchrotron radiation 
from the infalling gas in the same way as done previously 
with affine parameter. To create an image, one specifies 
a grid of points in a, (3 and calculates and I. This fully 
specifies the geodesic, and we can calculate the spacetime 
coordinates at which it intersects the accreting gas. The 
intensity along each geodesic represents a point in the 
image, which is why it is so important to be able to 
calculate geodesies rapidly. 

The redshift is calculated using Viergutz (1993). Here, 
the flow is spherically symmetric and. 
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with 



e-'^ = Ap^E-\ e2'*i=p2A-i, 7 = (1 - t;'-2)-i/2, 

(85) 

and v'^ is the radial component of the four- velocity. Writ- 
ten in terms of u in the Schwarzschild metric, this sim- 
plifies to 



VI -2u 



(86) 



When the w-component of the four- velocity, u", vanishes, 
this reproduces the standard gravitational redshift (Har- 
tlc 2003). 

We first ignore absorption and compare radiated spec- 
tra with the analytic calculation. The result is in Fig. 8. 
Shapiro (1973b) points out that the synchrotron radia- 
tion is dominated by a thin spherical shell of gas with 

~ i^c- Then the first part of the spectrum, where 
L^g ~ j/^/"^, originates from the outer part of the sphere. 
The bending of light should be negligible in that region 
and the ray tracing should agree with the analytic result, 
which it does to within ~ 5%. At higher frequencies, the 
radiation is originating in the innermost radii, and the 
bending of light becomes significant. The difference is 
~ 15% at high frequencies. 

Next, absorption is included. Fig. 9 compares the spec- 
tra calculated with and without absorption. The radia- 
tion is heavily attenuated at frequencies < 10^^ Hz. At 
these frequencies, the luminosity is dominated by the in- 
nermost optically thin radius, which we take to be the 
radius where r = 1. Blackbody emission at the tem- 
perature of gas at this radius, converted to a luminosity 
by integrating over impact parameter, is labeled 'Ther- 
mal' in Fig. 9 and is a decent approximation to the full 
spectrum when the fiuid is optically thick. 

Prom i/Q — 10® Hz to 1^0 — 10 Hz, the gas is op- 
tically thick everywhere. Then only thermal emission 
from the outermost radius is seen, and the spectrum 
follows a Rayleigh- Jeans curve with L^.^ ^ Uq. Prom 
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vq ~ 10^" Hz to — 10^^, the innermost optically thin 
radius is changing, and the luminosity begins to turn 
over. Starting at Vq ~ 10^^ Hz, the gas is optically thin to 
the synchrotron radiation, and the spectrum reduces to 
that of emission only (labeled 'Emission' in Fig. 9). This 
result agrees reasonably well with the assertion made 
by Shapiro (1973b) that absorption is negligible when 
u > 10" Hz. 

Also of interest is the black hole shadow produced 
by various accretion models (Falcke et al. 2000). Fig. 
10 shows the shadow of the spherically accreting 

Schwarzschild black hole as a 2-d contour plot and a 1-d 
profile. The shadow is produced at + ~ 27, and 
is caused by the difTcrcncc in proper length of geodesies 
which intersect the horizon and return to infinity, as well 
as the blueshift of radiation from infalling gas behind the 
black hole relative to the red shift of that nearest the ob- 
server. The asymmetry in Fig. 2 is not seen here due to 
the spherical symmetry of the Schwarzschild metric. 

8. FUTURE WORK 

The code presented here is the first to calculate all co- 
ordinates of Kerr null geodesies semi-analytically. This 
work's natural extension is to timelike geodesies, which 
involves many more cases, but only straightforward gen- 
eralizations of the formulas given here (RB94, Appendix 
A) . The main challenge is that for bound orbits it is dif- 
ficult to specify the number of radial turning points in 
advance. Ideally the afiine parameter could be used as 
an independent variable to indicate how far along the 
geodesic to trace. However, it is a function of u and /x 
which cannot be inverted. 

8.1. Advantages of Analyticity 

The main advantages of using a semi-analytic code 
such as that presented here for tracing geodesies are 
speed, accuracy and flexibility. The speed increase from 
our code depends greatly on the application considered. 
For ray tracing applications, a lower bound is a factor of 
5 in the case where all coordinates are being calculated, 
and the entire ray is being traced. The maximum speed 
increase is probably a factor between 100-500 in the case 
where the code is solving for geodesic coordinates at a 
specific point. 

The importance of speed in tracing geodesies depends 
on the computational expense of their construction rela- 
tive to that of the rest of the desired calculation. For the 
simple radiative transfer applications considered here, 
time spent computing geodesies dominates in creating 
Figs. 2,3,4. The construction of geodesies and radiative 
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transfer parts are about equally expensive in creating 
Figs. 10 and geodesic speed is relatively tmimportant in 
the calculations leading to Figs. 5-9. In the latter cases, 
this is because the same geodesies can be re-used at many 
time steps, frequencies or both. 

The trend from these toy problems is that the sim- 
pler cases benefit most from rapid geodesic calculation. 
However, there is reason to expect that for more realistic 
calculations rapid geodesic construction will again be im- 
portant. Most accretion flows transition from optically 
thin to thick. To accurately compute radiative transfer 
from such flows, it is often necessary to take small steps 
in the vicinity where the optical depth is about unity. 
This requires calculating extra geodesic trajectories in 
this region. Since the region where the optical depth 
changes rapidly depends on frequency, and for a time- 
dependent accretion model on time as well, new geodesies 
must be computed at each time step and frequency, and 
hence they cannot be re-used as is the ease for the time- 
independent, optically thin models considered in almost 
all our examples. 

The precision of our code is also extremely high over 
a broad range of geodesic parameters. This is currently 
less important in radiative transfer applications where 
the dynamical models are uncertain, but it is important 
in caustic calculations such as those in RB94 and Bozza 
(2008). Finally, our code can compute arbitrary sections 
of geodesies in any direction. This flexibility allows extra 
points to be calculated in regions where the optical depth 
is changing rapidly or to check convergence on the fly. It 
may also be useful in a future method for computing 
Compton scattering, in which rays are traced outwards 
from each point on the geodesic to calculate the scattered 
intensity into that point. 

Unlike previous analytic work, our code makes no as- 
sumption of time-independence or axisymmetry in the 
accretion flow and is therefore well suited to the geome- 
tries used in 3D GRMHD simulations. Computationally 
expensive observables such as polarization and variability 
will be much more tractable given the speed and flexibil- 
ity of this code. 
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